R Notebook

source('../settings/settings.R')
source('commonFunctions.R')

persons <- SELECTED_SUBJECTS
drive <- 4

inputFile <- str_interp('../data/processed/distancewise/TT1_Drive_${drive}_${distPrev}m_${distNext}m.csv', list(drive=drive, distPrev=30, distNext=30))
outputFile <- str_interp("../data/processed/analysis/TT1_Drive_${drive}_PP_${distPrev}m_${distNext}m.csv", list(drive=drive, distPrev=30, distNext=30))

all_Drive4 <- read.csv(inputFile)
all_Drive4$Subject <- as.factor(all_Drive4$Subject)
all_Drive4$logPerspiration <- log(all_Drive4$Perspiration)

# starting_points = c( 669 , 668 , 676 , 687 , 680 ,  676 ,  678 ,
#                      693 , 722 , 723 , 677 , 679 ,  711 ,  707 ,  
#                      699 , 679 , 684 , 688 , 686 ,  696 ,  702 )
# 
# ending_points   = c( 741 , 786 , 749 , 782 , 736 ,  756 ,  768 ,  
#                      812 , 853 , 792 , 783 , 772 ,  799 ,  781 ,  
#                      777 , 763 , 795 , 791 , 832 ,  755 ,  758 )

peak_points     = c( 67 ,  86 ,  73 ,  73 ,  73 ,  64 ,  73 ,  
                     79 ,  69 ,  64 ,  68 ,  67 ,  77 ,  68 ,  
                     82 ,  67 ,  72 ,  72 ,  71 ,  68 ,  64 )

# Driving time
driving_times = vector(mode="list", length = length(persons))
names(driving_times) <- persons

activity_names = vector(mode="list", length = length(persons))
names(activity_names) <- persons

acc_start_times = vector(mode="list", length = length(persons))
names(acc_start_times) <- persons
acc_end_times = vector(mode="list", length = length(persons))
names(acc_end_times) <- persons

stressor_start_times = vector(mode="list", length = length(persons))
names(stressor_start_times) <- persons
stressor_end_times = vector(mode="list", length = length(persons))
names(stressor_end_times) <- persons

complete_times = vector(mode="list", length = length(persons))
names(complete_times) <- persons

data_baseline = vector(mode="list", length=length(persons))
pp_baseline = vector(mode="list", length=length(persons))

names(data_baseline) <- persons
names(pp_baseline) <- persons

# Number of peaks

PREV_DISTANCE = 300
TRACKING_DISTANCE = 60
DRIVE_MODE = 4

getActivityName <- function(x, fullname=F) {
  if(x == 1) return(ifelse(fullname, "Normal", "NO"))
  if(x == 2) return(ifelse(fullname, "Cognitive", "C"))
  if(x == 3) return(ifelse(fullname, "Motoric", "M"))
}

for (p in persons) {
  pData <- all_Drive4[all_Drive4$Subject==as.integer(p) | all_Drive4$Subject==p,]
  pAcc <- pData[pData$Failure>0.5,] # Failure = 1
  acc_start_times[[p]] <- min(pAcc$Distance)
  acc_end_times[[p]] <- max(pAcc$Distance)
  
  activity_names[[p]] <- getActivityName(pData[pData$Time==60,]$Activity, fullname = T)
  
  pStressor = pData[pData$Activity>1.5,] # Stressor = 2, 3
  if (nrow(pStressor) > 0) {
    stressor_start_times[[p]] <- min(pStressor$Distance)
    stressor_end_times[[p]] <- max(pStressor$Distance)
  } else {
    stressor_start_times[[p]] <- NULL
    stressor_end_times[[p]] <- NULL
  }
}

DELAY_DISTANCE <- 10

idx <- 1
plt_AllAcc <- vector(mode="list", length=length(persons)) 
names(plt_AllAcc) <- persons

COLOR_ACC = "#02A3C8"
COLOR_PP = "#F28E8E"
COLOR_BRAKE = "#888888"

y1 <- list(
  tickfont = list(color = COLOR_ACC),
  title="Degree",
  range=c(0, 100)
)
y2 <- list(
  tickfont = list(color = COLOR_PP),
  overlaying = "y",
  side = "right",
  title = "Log Perspiration",
  showgrid = FALSE,
  range=c(min(all_Drive4$ppLogNormalized), max(all_Drive4$ppLogNormalized))
)
  
for (p in persons) {
  pData <- all_Drive4[all_Drive4$Subject==as.integer(p) | all_Drive4$Subject==p,]
  
  # Baseline
  data_baseline[[p]] <- read.csv(str_interp("../data/processed/drives/T0${person}/T0${person}_Drive_1.csv", list(person=p)))
  # Compute the mean
  p_pp_nr <- data_baseline[[p]]$Perspiration
  p_pp_nr <- p_pp_nr[!is.na(p_pp_nr)]
  pp_baseline[[p]] <- log(mean(p_pp_nr))
  
  # Incident
  driving_times[[p]] <- max(pData$Distance)
  
  incident_starting_time <- acc_start_times[[p]] # starting_points[idx]
  incident_ending_time <- acc_end_times[[p]]
  complete_times[[p]] <- ifelse(incident_starting_time + TRACKING_DISTANCE > driving_times[[p]], driving_times[[p]], incident_starting_time + TRACKING_DISTANCE)
  
  from_time <- ifelse(incident_starting_time - PREV_DISTANCE >= 0, incident_starting_time - PREV_DISTANCE, 0)
  to_time <- complete_times[[p]]
  
  # print(paste("From", from_time))
  # print(paste("Incident", incident_starting_time))
  # print(paste("To", to_time))
  
    
  pDataBefore <- pData[pData$Distance < incident_starting_time & pData$Distance >= from_time,]
  pDataAfter <- pData[pData$Distance >= incident_starting_time + DELAY_DISTANCE & pData$Distance <= to_time,]
  
  # print(nrow(pDataBefore))
  # print(nrow(pDataAfter))
  
  ppMeanBefore <- mean(pDataBefore$ppLogNormalized)
  ppMeanAfter <- mean(pDataAfter$ppLogNormalized)
  
  dir.create(file.path('../figures/drive/', paste0('Drive_', DRIVE_MODE)), showWarnings = FALSE)
  fname <- str_interp('../figures/drive/Drive_${drive}/P${person}.svg', list(drive=DRIVE_MODE, person=p)) 
  
  pData <- pData[pData$Distance >= from_time,]
  plot_Acc <- plot_ly(pData, x = ~Distance, height=400, width=900) %>%
              add_trace(name="Acceleration", y = ~Acceleration, type = 'scatter', mode = 'lines', line=list(width=1.5, color=COLOR_ACC)) %>% 
              add_trace(name="Brake", y = ~Braking, type = 'scatter', mode = 'lines', line=list(width=1.5, color=COLOR_BRAKE)) %>%
              add_trace(name="PP", y = ~ppLogNormalized, type = 'scatter', mode = 'lines', line=list(width=1.5, color=COLOR_PP), yaxis = "y2") %>% 
              add_segments(x = min(pData$Distance), xend = max(pData$Distance), y = ppMeanBefore, yend = ppMeanBefore, 
                           yaxis = "y2", name="Mean PP (Before Incident)",
                           line=list(color=COLOR_PP, dash = 'dot')) %>%
              add_segments(x = min(pData$Distance), xend = max(pData$Distance), y = ppMeanAfter, yend = ppMeanAfter, 
                           yaxis = "y2", name="Mean PP (After Incident)",
                           line=list(color="darkred", dash = 'dot')) %>%
              # add_segments(x = min(pData$Distance) - 0.1, xend = max(pData$Distance), y = pp_baseline[[p]], yend = pp_baseline[[p]], 
              #              yaxis = "y2", name="Baseline PP (from Drive 1)",
              #              line=list(color="blue", dash = 'dot')) %>%
    
              layout(
                title=paste0("Subject #", p, " (Stressor=", activity_names[[p]], ")"), 
                xaxis=list(title="Distance [m]", range=c(0)), 
                yaxis=y1, 
                yaxis2=y2, 
                margin = list(l = 50, r = 50, b = 50, t = 50, pad = 4),
                shapes = list(
                  # Holistic period
                  list(type = "rect", fillcolor = "red", 
                       line = list(color = "red"), opacity = 0.3,
                      x0 = incident_starting_time, x1 = incident_ending_time, xref = "x",
                      y0 = 0, y1 = 100, yref = "y"),
                  # Stressor period
                  list(type = "rect", fillcolor = "yellow", 
                       line = list(color = "yellow"), opacity = 0.1,
                      x0 = stressor_start_times[[p]], x1 = incident_starting_time, xref = "x",
                      y0 = 0, y1 = 100, yref = "y"),
                  list(type = "rect", fillcolor = "yellow", 
                       line = list(color = "yellow"), opacity = 0.1,
                      x0 = incident_ending_time, x1 = stressor_end_times[[p]], xref = "x",
                      y0 = 0, y1 = 100, yref = "y")
                ),
                legend = list(x = 0.1, y = 1, bgcolor = "rgba(0,0,0,0)", title="Metric"),
                autosize = F
              )
  
  # orca(plot_PP, fname)
  idx <- idx + 1
  plt_AllAcc[[p]] <- plot_Acc
}

htmltools::tagList(plt_AllAcc)

# for (p in persons) {
#   # Save image
#   orca(plt_AllAcc[[p]], file = paste0("../plots/drive/Drive_4/T0", p, ".png"), scale = 2)
# }

idx <- 1
behavioralColumns <- c("Subject", 
                       "Brake_u", 
                       "Brake_std", 
                       "PP_before",
                       "PP_u",  
                       "PP_std",
                       "PP_dev")
behavioralMatrix <- matrix(nrow=length(persons), ncol = length(behavioralColumns))

# Careful about Subject 09
# selected_persons <- persons[persons != "09"]

for (p in persons) {
  pData <- all_Drive4[all_Drive4$Subject==as.integer(p) | all_Drive4$Subject==p,]
  
  incident_starting_time <- acc_start_times[[p]] # starting_points[idx]
  incident_ending_time <- acc_end_times[[p]]
  
  from_time <- ifelse(incident_starting_time - PREV_DISTANCE >= 0, incident_starting_time - PREV_DISTANCE, 0)
  to_time <- complete_times[[p]]
    
  dfBefore <- pData[pData$Distance < incident_starting_time & pData$Distance >= from_time,]
  dfAfter <- pData[pData$Distance >= incident_starting_time + DELAY_DISTANCE  & pData$Distance <= to_time,]
  
  # diffSpeed <- mean(dfAfter$Speed) - mean(dfBefore$Speed)
  brakeMean <- mean(dfAfter$Braking)
  brakeStd <- mean(dfAfter$Braking)
  
  ppMean <- mean(dfAfter$ppLogNormalized)
  ppBefore <- mean(dfBefore$ppLogNormalized)
  ppStd <- sd(dfAfter$ppLogNormalized)
  
  mid_avg <- (pp_baseline[[p]] + mean(dfBefore$ppLogNormalized)) / 2
    
  diffPP <- mean(dfAfter$ppLogNormalized) - mean(dfBefore$ppLogNormalized)
  
  behavioralMatrix[idx, ] <- c(p, 
                               round(brakeMean, digits=5), 
                               round(brakeStd, digits=5),
                               round(ppBefore, digits=5),
                               round(ppMean, digits = 5),
                               round(ppStd, digits=5),
                               round(diffPP, digits=5))
  idx <- idx + 1
}

# behavioralMatrix

behavioralDf <- as.data.frame(behavioralMatrix, stringsAsFactors=FALSE)
names(behavioralDf) <- behavioralColumns

behavioralDf

NA

clusteringDf <- behavioralDf
clusteringDf$Subject <- NULL
# clusteringDf$PP_dev_norm <- as.numeric(clusteringDf$PP_dev) / as.numeric(clusteringDf$PP_u)
# clusteringDf$PP_std_norm <- as.numeric(clusteringDf$PP_std) / as.numeric(clusteringDf$PP_u)
clusteringDf$Brake_u <- NULL
clusteringDf$Brake_std <- NULL
clusteringDf$PP_before <- NULL
clusteringDf$PP_u <- NULL
clusteringDf$PP_std <- NULL
# clusteringDf$PP_dev <- NULL

rownames(clusteringDf) <- paste0("#", persons)

for (col in names(clusteringDf)) {
  clusteringDf[,col] <- as.numeric(as.character(clusteringDf[, col]))
  # clusteringDf[,col] <- scale(clusteringDf[,col])
}
clusteringDf

dfActivity <- all_Drive4[all_Drive4$Time==60,] %>% select(c("Subject", "Activity"))
dfActivity$ActivityName <- sapply(dfActivity$Activity, getActivityName)
dfActivity$Subject <- as.factor(dfActivity$Subject)
rownames(dfActivity) <- NULL
dfActivity %>% select(c("Subject", "ActivityName"))

library(dendextend)

NUMBER_OF_CLUSTERS = 3

color_darkpink = "#e75480"
CLUSTER_BRANCH_COLORS <- c("blue", "red", color_darkpink)[1:NUMBER_OF_CLUSTERS]
CLUSTER_LABEL_COLORS <- c("blue", "red", color_darkpink)[1:NUMBER_OF_CLUSTERS]

behavioralMatrixClustering <- as.matrix(clusteringDf)
rownames(behavioralMatrixClustering) <- paste0(dfActivity$ActivityName, " - #", persons)
distMatrix <- dist(behavioralMatrixClustering, method="manhattan")
hresults <- distMatrix %>% hclust

hc <- hresults %>% 
      as.dendrogram %>%
      set("nodes_cex", NUMBER_OF_CLUSTERS) %>%
      set("labels_col", value = CLUSTER_LABEL_COLORS, k=NUMBER_OF_CLUSTERS) %>%
      # set("leaves_pch", 19) %>%
      # set("leaves_col", value = c("gray"), k=NUMBER_OF_CLUSTERS) %>%    
      set("branches_k_color", value=CLUSTER_BRANCH_COLORS, k=NUMBER_OF_CLUSTERS)

plot(hc)
legend("topright", 
     title="Drive=Failure \nChange of Arousal",
     legend = c("Exceptional Increase" , "Noticable Increase" , "No-change or Decrease"), 
     col = c("red", "pink" , "blue"),
     pch = c(20,20,20), bty = "n",  pt.cex = 1.5, cex = 0.8 , 
     text.col = "black", horiz = FALSE, inset = c(0.0, 0.1))

NUMBER_OF_CLUSTERS = 2

color_darkpink = "#e75480"
CLUSTER_BRANCH_COLORS <- c("blue", "red", color_darkpink)[1:NUMBER_OF_CLUSTERS]
CLUSTER_LABEL_COLORS <- c("blue", "red", color_darkpink)[1:NUMBER_OF_CLUSTERS]

behavioralMatrixClustering <- as.matrix(clusteringDf)
rownames(behavioralMatrixClustering) <- paste0(dfActivity$ActivityName, " - #", persons)
distMatrix <- dist(behavioralMatrixClustering, method="manhattan", diag = T)
hresults <- distMatrix %>% hclust(method="average")

hc <- hresults %>% 
      as.dendrogram %>%
      set("nodes_cex", NUMBER_OF_CLUSTERS) %>%
      set("labels_col", value = CLUSTER_LABEL_COLORS, k=NUMBER_OF_CLUSTERS) %>%
      # set("leaves_pch", 19) %>%
      # set("leaves_col", value = c("gray"), k=NUMBER_OF_CLUSTERS) %>%    
      set("branches_k_color", value=CLUSTER_BRANCH_COLORS, k=NUMBER_OF_CLUSTERS)

plot(hc)
legend("topright", 
     title="Drive=Failure \nChange of Arousal",
     legend = c("Accelerophobia" , "Normal"), 
     col = c("red", "blue"),
     pch = c(20,20,20), bty = "n",  pt.cex = 1.5, cex = 0.8 , 
     text.col = "black", horiz = FALSE, inset = c(0.0, 0.1))

# Store clustering data
dfx <- clusteringDf
dfx <- cbind(persons, as.numeric(behavioralDf$PP_before), as.numeric(behavioralDf$PP_u), dfx, dfActivity$ActivityName)
names(dfx) <- c("Subject", "PP_Prior", "PP_After", "PP_Dev", "Activity")

# Normalize
dfx <- dfx %>% mutate(PP_Dev_Normalized = ifelse(PP_Prior > 0, PP_After, PP_After - PP_Prior))
write.csv(dfx, outputFile, row.names = F)

library(cluster)
fit <- kmeans(clusteringDf, 2)
clusplot(clusteringDf, fit$cluster, color=TRUE, shade=TRUE,
   labels=2, lines=0)

silhouette_score <- function(k){
  km <- kmeans(clusteringDf, centers = k, nstart=25)
  ss <- silhouette(km$cluster, dist(clusteringDf))
  mean(ss[, 3])
}
k <- 2:10
avg_sil <- sapply(k, silhouette_score)
plot(k, type='b', avg_sil, xlab='Number of clusters', ylab='Average Silhouette Scores', frame=FALSE)

Linear Model

sampledData <- getSampleSegmentedData(NA, all_Drive4, window=2)
linearModelOnline <- lmer(ppNext ~ 
              (1 | Subject)
              + Speed_u
              + Speed_std
              + Acc_u
              + Acc_std
              + Brake_u
              + Brake_std
              + Steering_u
              + Steering_std, 
            data=sampledData, REML = T)

# lmer(ppLogNormalized ~ (1 | Subject) + Speed_u + Speed_std + Acc_u + Acc_std + Brake_u + Brake_std + Steering_u + Steering_std + HR + BR, data = pData, REML = T)

# anova(model)
summary(linearModelOnline)

Linear mixed model fit by REML. t-tests use Satterthwaite's method ['lmerModLmerTest']
Formula: ppNext ~ (1 | Subject) + Speed_u + Speed_std + Acc_u + Acc_std +      Brake_u + Brake_std + Steering_u + Steering_std
   Data: sampledData

REML criterion at convergence: -1352.9

Scaled residuals: 
    Min      1Q  Median      3Q     Max 
-2.4668 -0.6034 -0.0884  0.5493  5.5490 

Random effects:
 Groups   Name        Variance Std.Dev.
 Subject  (Intercept) 0.015927 0.12620 
 Residual             0.007908 0.08893 
Number of obs: 775, groups:  Subject, 21

Fixed effects:
               Estimate Std. Error         df t value Pr(>|t|)    
(Intercept)   3.994e-02  3.124e-02  3.184e+01   1.278   0.2103    
Speed_u      -5.214e-04  5.083e-04  7.471e+02  -1.026   0.3053    
Speed_std     5.737e-03  2.795e-03  7.467e+02   2.053   0.0405 *  
Acc_u        -1.140e-04  4.334e-04  7.472e+02  -0.263   0.7925    
Acc_std       1.666e-03  8.636e-04  7.463e+02   1.929   0.0541 .  
Brake_u       3.263e-03  4.323e-04  7.469e+02   7.548 1.29e-13 ***
Brake_std     3.042e-04  7.076e-04  7.482e+02   0.430   0.6673    
Steering_u   -2.200e-04  5.226e-05  7.506e+02  -4.211 2.85e-05 ***
Steering_std -9.907e-05  1.384e-04  7.472e+02  -0.716   0.4744    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Correlation of Fixed Effects:
            (Intr) Speed_ Spd_st Acc_u  Acc_st Brake_ Brk_st Sterng_
Speed_u     -0.409                                                  
Speed_std   -0.027  0.054                                           
Acc_u       -0.129 -0.155 -0.147                                    
Acc_std      0.081 -0.043 -0.052 -0.494                             
Brake_u     -0.131  0.145 -0.657  0.316 -0.077                      
Brake_std   -0.275  0.563  0.328  0.032 -0.149 -0.312               
Steering_u  -0.008  0.009  0.095 -0.013  0.019 -0.079 -0.036        
Steerng_std -0.100  0.206 -0.398  0.110 -0.204  0.385 -0.349 -0.058 

plot(linearModelOnline)